From a 1m DSM in Tiraumea, New Zealand, terrain derivatives will be computed with the following code. The aim is to detect Earthflows, possibly through machine learning, where a model would be fed with the terrain derivatives, and optical information.
Since the data is large, two initial study areas were proposed for analysis:
# Call mapping library and set options
library(mapview)
mapviewOptions(basemaps = 'OpenStreetMap', console = F, verbose = F, leafletWidth = 800)
library(sf, quietly = T, verbose = F, warn.conflicts = F)
tiraumea = st_read('study_area/Tiraumea.shp', quiet = T)
sa = st_read('study_area/study_area_outlines.shp', quiet = T)
mapview(list(tiraumea,sa))
Initially, we will look at Study Area 1. I will load the DEM data into R with the help of the raster package. The DEM sinks were previously filled with this example code:
rsaga.fill.sinks(in.dem = "dem.sgrd", out.dem = "dem_filled", method = "wang.liu.2006", env = env)
Then it was cropped to the study area (in QGIS).
Now for the resulting DEM:
# I use 'raster` at this point, since for some reason writing a `stars` object into a SAGA compatible format is not working.
library(raster, quietly = T)
dem = raster('input_dem/dsm_filled_sa1.tif')
crs(dem) = '+init=epsg:2193'
mapview(dem)
I save it as a SAGA compatible format to be able to call it into the RSAGA functions.
writeRaster(dem, 'input_dem/dsm_sa1.sgrd', format = 'SAGA', overwrite = T, NAflag = 0, prj = T)
First I set up the RSAGA environment:
library(RSAGA)
#library(link2GI)
#saga = linkSAGA()
#saga # Copy the results of this to the env setting
env = rsaga.env(
path = "C:\\Users\\b1066081\\Desktop\\saga-7.5.0_x64",
modules = "C:\\Users\\b1066081\\Desktop\\saga-7.5.0_x64\\tools"
# cmd = "C:\\Users\\b1066081\\Desktop\\saga-7.5.0_x64\\saga_cmd.exe"
)
Verify specified path to SAGA command line program...
Verify specified path to SAGA modules...
Done
And now, compute some derivatives for the different modules.
With this I will compute overall 8 outputs:
I will also compute variables to keep with the workflow of Eisank et al. (2014), but I will not add the multiscale parameter, since I do not find how this was done for that.
rsaga.slope.asp.curv(
in.dem = "input_dem/dsm_sa1.sgrd",
out.slope = "out_products/slope_sa1.sgrd",
out.cgene = "out_products/cgene_sa1.sgrd",
out.cplan = "out_products/cplan_sa1.sgrd",
out.cprof = "out_products/cprof_sa1.sgrd",
out.ccros = "out_products/ccros_sa1.sgrd",
out.clong = "out_products/clong_sa1.sgrd",
out.cmaxi = "out_products/cmaxi_sa1.sgrd",
out.cmini = "out_products/cmini_sa1.sgrd",
unit.slope = "degrees", env = env,
method = "poly2zevenbergen"
)
From this modules on, I will need to find out the usage first, because there are no shortcuts for them. I will always refer to this table for the morphometry module to get the respective code:
rsaga.get.modules('ta_morphometry', env = env) %>% knitr::kable()
|
rsaga.get.usage('ta_morphometry', 14, env = env)
From here I will get 5 outputs, and will work with default W, T, E.
rsaga.geoprocessor(
'ta_morphometry', 14,
list(
DEM='input_dem/dsm_sa1.sgrd',
HO='out_products/slhgt_sa1.sgrd',
HU='out_products/vldpt_sa1.sgrd',
NH='out_products/nrhgt_sa1.sgrd',
SH='out_products/sthgt_sa1.sgrd',
MS='out_products/mdslp_sa1.sgrd'
),
env = env
)
This tool has its own module, and will return one outcome, coded ‘ddgrd’. We get the usage:
rsaga.get.usage('ta_morphometry', 9, env = env)
And call it, with default distance (10) and gradient in degrees:
rsaga.geoprocessor(
'ta_morphometry', 9,
list(
DEM='input_dem/dsm_sa1.sgrd',
GRADIENT='out_products/ddgrd_sa1.sgrd',
OUTPUT = 2
),
env = env
)
This tool has its own module, and will return one outcome, coded ‘mbidx’. We get the usage:
rsaga.get.usage('ta_morphometry', 10, env = env)
And compute with all the defaults:
rsaga.geoprocessor(
'ta_morphometry', 10,
list(
DEM='input_dem/dsm_sa1.sgrd',
MBI='out_products/mbidx_sa1.sgrd'
),
env = env
)
This tool has its own module, and will return one outcome, coded ‘tpidx’. We get the usage:
rsaga.get.usage('ta_morphometry', 18, env = env)
And compute with all the defaults, min and max are definitely needed and vary depending on the raster resolution. If not used properly computation will be exhaustive.
rsaga.geoprocessor(
'ta_morphometry', 18,
list(
DEM='input_dem/dsm_sa1.sgrd',
TPI='out_products/tpidx_sa1.sgrd',
RADIUS_MIN=0,
RADIUS_MAX=20
),
env = env
)
This tool has its own module, and will return one outcome, coded ‘cvidx’. We get the usage:
rsaga.get.usage('ta_morphometry', 1, env = env)
And compute with method gradient and neighbours 3x3:
rsaga.geoprocessor(
'ta_morphometry', 1,
list(
ELEVATION='input_dem/dsm_sa1.sgrd',
RESULT='out_products/cvidx_sa1.sgrd',
METHOD=1,
NEIGHBOURS=1
),
env = env
)
This tool has its own module, and will return one outcome, coded ‘tridx’. We get the usage:
rsaga.get.usage('ta_morphometry', 16, env = env)
And compute with all the defaults:
rsaga.geoprocessor(
'ta_morphometry', 16,
list(
DEM='input_dem/dsm_sa1.sgrd',
TRI='out_products/tridx_sa1.sgrd'
),
env = env
)
This tool has its own module, and will return one outcome, coded ‘textu’. We get the usage:
rsaga.get.usage('ta_morphometry', 20, env = env)
And compute with all the defaults:
rsaga.geoprocessor(
'ta_morphometry', 20,
list(
DEM='input_dem/dsm_sa1.sgrd',
TEXTURE='out_products/textu_sa1.sgrd'
),
env = env
)
Again I will access the codes for the tools to get the usage:
rsaga.get.modules('ta_hydrology', env = env) %>% knitr::kable()
|
This tool has its own module, and will return one outcome, coded ‘mridx’. We get the usage:
rsaga.get.usage('ta_hydrology', 23, env = env)
And compute with all the defaults:
rsaga.geoprocessor(
'ta_hydrology', 23,
list(
DEM='input_dem/dsm_sa1.sgrd',
MRN='out_products/mridx_sa1.sgrd'
),
env = env
)
This particular tool I had to compute using the GUI for SAGA, because RSAGA does not let me use the LS Factor (one step) library that I need to do it right. In the GUI I only need to give the DEM. It is saved as ‘lsfct’.
This tool has its own module, and will return one outcome, coded ‘sagaw’. We get the usage:
rsaga.get.usage('ta_hydrology', 15, env = env)
And compute with defaults:
rsaga.geoprocessor(
'ta_hydrology', 15,
list(
DEM='input_dem/dsm_sa1.sgrd',
TWI='out_products/sagaw_sa1.sgrd'
),
env = env
)
This tool has its own module, and will return one outcome, coded ‘sllgt’. We get the usage:
rsaga.get.usage('ta_hydrology', 7, env = env)
And compute with defaults:
rsaga.geoprocessor(
'ta_hydrology', 7,
list(
DEM='input_dem/dsm_sa1.sgrd',
LENGTH='out_products/sllgt_sa1.sgrd'
),
env = env
)
To compute this index, some previous steps are required. The version incompatibility of SAGA 7.5 and RSAGA does not support this computation(?), at least in an intuitive way. So I compute on the SAGA GUI.
The steps I took were:
ta_hidrology module. The only input is the DEM. The intermediate output needed is the Specific Catchment Area. I saved this on the int_products directory as ‘spcar_sa1’. Specific parameters are:
ta_hidrology module. Inputs are the Slope and Specific Catchment Area. The result is saved to out_products as ‘spidx’.This particular tool I had to compute using the GUI for SAGA, because RSAGA does not let me use the TWI (one step) library that I need to do it right. In the GUI I only need to give the DEM. I selected Deterministic 8 for the Flow Distribution. It is coded ‘twidx’.
Again I will access the codes for the tools to get the usage:
rsaga.get.modules('ta_lighting', env = env) %>% knitr::kable()
|
For this and the next tool, a multiscale parameter is chosen to speed up the computation. From this tool we can get two derivatives:
rsaga.get.usage('ta_lighting', 3, env = env)
And compute with defaults:
rsaga.geoprocessor(
'ta_lighting', 3,
list(
DEM='input_dem/dsm_sa1.sgrd',
VISIBLE='out_products/visky_sa1.sgrd',
SVF='out_products/svfct_sa1.sgrd',
METHOD=1,
DLEVEL=3
),
env = env
)
From this tool we can get two derivatives:
rsaga.get.usage('ta_lighting', 5, env = env)
library path: C:\Users\b1066081\Desktop\SAGA-7~1.0_X\tools\
library name: ta_lighting
library : ta_lighting
Usage: saga_cmd ta_lighting 5 [-DEM <str>] [-POS <str>] [-NEG <str>] [-RADIUS <double>] [-METHOD <str>] [-DLEVEL <double>] [-NDIRS <num>] [-UNIT <str>] [-NADIR <str>]
-DEM:<str> Elevation
Grid (input)
-POS:<str> Positive Openness
Grid (output)
-NEG:<str> Negative Openness
Grid (output)
-RADIUS:<double> Radial Limit
Floating point
Minimum: 0.000000
Default: 10000.000000
-METHOD:<str> Method
Choice
Available Choices:
[0] multi scale
[1] line tracing
Default: 1
-DLEVEL:<double> Multi Scale Factor
Floating point
Minimum: 1.250000
Default: 3.000000
-NDIRS:<num> Number of Sectors
Integer
Minimum: 2
Default: 8
-UNIT:<str> Unit
Choice
Available Choices:
[0] Radians
[1] Degree
Default: 0
-NADIR:<str> Difference from Nadir
Boolean
And compute with defaults:
rsaga.geoprocessor(
'ta_lighting', 5,
list(
DEM='input_dem/dsm_sa1.sgrd',
POS='out_products/posop_sa1.sgrd',
NEG='out_products/negop_sa1.sgrd',
METHOD=0,
DLEVEL=3
),
env = env
)
For this some previous steps are needed, see the complete code:
rsaga.get.usage('ta_hydrology', 0, env = env)
rsaga.geoprocessor(
'ta_hydrology', 0,
list(
ELEVATION='input_dem/dsm_sa1.sgrd',
FLOW='int_products/flowsa1.sgrd',
METHOD=0
),
env = env
)
rsaga.get.usage('ta_channels', 0, env = env)
rsaga.geoprocessor(
'ta_channels', 0,
list(
ELEVATION='input_dem/dsm_sa1.sgrd',
CHNLNTWRK='int_products/chnet_sa1.sgrd',
INIT_GRID='int_products/flowsa1.sgrd',
INIT_VALUE=1000000
),
env = env
)
rsaga.get.usage('ta_channels', 3, env = env)
rsaga.geoprocessor(
'ta_channels', 3,
list(
ELEVATION='input_dem/dsm_sa1.sgrd',
CHANNELS='int_products/chnet_sa1.sgrd',
DISTANCE='out_products/vdcnw_sa1.sgrd'
),
env = env
)
Useful function to convert to .tif extension.
files = list.files(path = "./out_products", pattern = "sa1*.*sdat", full.names = T)
convertToTiff <- function(x) {
b = brick(x, crs = '+init=epsg:2193')
b = projectRaster(b, crs = '+init=epsg:2193')
extension(x) = '.tif'
dt = dataType(b)
writeRaster(b, filename = x, datatype = dt)
}
lapply(files, convertToTiff)
And fast, view the results:
files = list.files(path = "./out_products", pattern = "*sa1.*tif", full.names = T)
file_names = list.files(path = "./out_products", pattern = "*sa1.*tif")
library(stars)
library(purrr)
files_read = files %>% map(function(x) read_stars(x, proxy = T))
par(mfrow = c(5,6))
invisible(lapply(files_read, function(x) {
plot(x, key.pos = NULL, reset = F, main = NULL)
}))
list = read.csv('code_list.csv', sep = ';')
list
Eisank, C., Hölbling, D., & Friedl, B. (2014). How well do terrain objects derived from pre-event DEM spatially correspond to landslides? Geological Society of America Abstracts with Programs, 46(6), 105. Vancouver, Canada.